	PROGRAM LACOM

C Gas phase chemistry by Liisa Pirjola (box model)
c for the Board of Education
c solar intensity is an array of 79*57 called aflux3.dat 
c zenith angle 30-86 astetta is enough in Finland

	IMPLICIT DOUBLE PRECISION ( A - Z )
	INTEGER N,MAX,I,IW,IFAIL,I1,II,I7,I4,KK,I33,par
	integer YEAR,day,mon,hour,pv,ku,tu,cl,nh,koodi,type1,nh1,lkm
	PARAMETER (N=67, IW=(12+N)*N+50)
	DOUBLE PRECISION TT,TEND,TOL
	DOUBLE PRECISION W(IW),Y(N)
	DOUBLE PRECISION J(79,57),ALFA
	COMMON /INTENS/ J,ALFA
	SAVE /INTENS/
	DOUBLE PRECISION K(123),KF(16),H2O,O2,N2,H2,M
	COMMON /LOPUT/ K,KF,H2O,O2,N2,H2,M
	SAVE /LOPUT/
	DOUBLE PRECISION V(67),EMI(67),VD(67),EMI1(67),wout
	COMMON /NOP/ V,EMI,VD,EMI1,wout
	SAVE /NOP/
	DOUBLE PRECISION T
	COMMON / TEMP /T
	SAVE / TEMP /
	EXTERNAL D02EAF,FCN

c simulation time (max), boundary layer height (h), latitude (ALAT), 
c longitude (ALON), starting day (DAY), month (MON), hour (HOUR), 
c albedo (ALB), uv radiation coefficient (G)
	OPEN(17,FILE='INPUT1.DAT', STATUS = 'OLD')
	READ(17,*) MAX,H,ALAT,ALON,YEAR,DAY,MON,HOUR,ALB,G,par
c ticu0 is the starting time from the beginning of the year
	ticu0=timecum(YEAR,day,mon,hour)
c        write(6,*) ticu0
	TEND=1.D-5
	TENDAPU=TEND
	TT = 0.
	RAJA=1.

C initial concentrations (molecule cm-3), deposition velocities (m s-1)
C  and emission rates (molecule cm-2 s-1)
      OPEN(16,FILE = 'INPUT2.DAT', STATUS='OLD')
      DO 15 I=1,54
       READ(16,*) I33,Y(I),VD(I),EMI1(I)
15    CONTINUE
      do 16 i=55,67
	y(i)=0
	vd(i)=0
	emi1(i)=0
16    continue

c calls subroutine which calculates the zenith angle
      CALL SOLAR(ALAT,ALON,YEAR,DAY,MON,HOUR)
      open(11,file='zeni.dat',status ='old')

C reads intensity of the radiation into an array (79x57), the intensity
c depends on the zenith angle and the wave length of the radiation
      OPEN(12,FILE = 'AFLUX3.DAT',STATUS = 'OLD')
      DO 300 KK=1,79
      READ(12,*) (J(KK,II), II=1,57)
300   CONTINUE
      CLOSE(12)
c the intensity is multiplied by albedo and uv-coefficient 
      if (alb.ge.0.46) then
	do 301 kk=1,79
	do 302 II=1,57
	J(KK,II)=(1+alb)*G*J(KK,II)
302   continue
301   continue
      else
       do 305 kk=11,79
       do 306 ii=1,57
       j(kk,ii)=(1+alb)*G*J(kk,ii)
306   continue
305   continue
      endif

c      write(6,*) '1+alb= ',1+alb,'uv coeff= ',g
C open the restore files and save the initial values
      OPEN(13,file='test1.res',status='unknown')
      OPEN(14,file='test2.res',status='unknown')
      WRITE(13,'(28e12.5)') TT,(Y(I1),I1=1,27)
      write(14,'(28e12.5)') TT,(Y(I1),I1=28,54)


C do 200 main-loop begins
      DO 200 I=1,MAX
      TICU=TICU0+TT/86400.
c      write(6,*) 'ticu = ',ticu
      OPEN(20,FILE='O3KRATE.DAT',STATUS='OLD')
      OPEN(21,FILE='NO2.DAT',STATUS='OLD')
      OPEN(22,FILE='NO3.DAT',STATUS='OLD')
      OPEN(24,FILE='HNO3.DAT',STATUS='OLD')
      OPEN(25,FILE='H2O2.DAT',STATUS='OLD')
      OPEN(26,FILE='HCHO.DAT',STATUS='OLD')
      OPEN(27,FILE='CH3CHO.DAT',STATUS='OLD')
      OPEN(28,FILE='CHOCHO.DAT',STATUS='OLD')
      OPEN(29,FILE='CH3OOH.DAT',STATUS='OLD')
      OPEN(30,FILE='CH3COCHO.DAT',STATUS='OLD')
	   
      read(11,*) pv,ku,tu,alfa
      alfa=real(nint(alfa))
c      write(6,*) 'alfa= ',alfa
      OPEN(10,FILE='TEMP.DAT', STATUS='OLD')
      OPEN(15,FILE='HUMI.DAT', STATUS='OLD')
      
c calculates the temperature in K
228   read(10,*) PV,KU,TU,TE2
c      write(6,*) 'min tai max l„mp”tila ',te2
      TICU2=TIMECUM(YEAR,PV,KU,TU)
      IF (TICU2.LE.TICU) THEN
	TICU1=TICU2
	TE1=TE2
	GOTO 228
      ENDIF
      close(10)
      T=TE1+(TE2-TE1)/12.*(TICU-TICU1)
      T=T+273.15
c      WRITE(6,*) 'l„mp”tila Kelvinein„= ',T

c calculates H2O-concentration in molecule cm-3
229   read(15,*) PV,KU,TU,RH1,TE1
      TICU1=TIMECUM(YEAR,PV,KU,TU)
      IF (TICU1.LE.TICU) THEN
c       write(6,*) 'ticu1 = ',ticu1
       RH=0.01*RH1
       TE=TE1+273.15
       GOTO 229
      ENDIF
c      WRITE(6,*) 'RH= ',RH,'l„mp”tila = ',TE
      psat = DEXP(77.344913-7235.4247/TE-8.2*DLOG(TE)+0.0057113*TE)
      H2O = psat*RH/(TE*1.381d-17)
c      WRITE(6,*) 'vesi = ',H2O
      close(15)

c calculates the cloudiness coefficient with which the intensity is 
c multiplied
	open(18,file='clou.dat', status='old')
230   read(18,*) pv,ku,tu,type1,nh1
	ticu1=timecum(YEAR,pv,ku,tu)
	if (ticu1.le.ticu) then
	  cl=type1
	  nh=nh1
	  goto 230
	endif
c      write(6,*) 'tyyppi= ',cl,'m„„r„= ',nh
      close(18)
      CALL CLOUDINESS(CL,NH,ALFA,C1,koodi)
      if (koodi.eq.1) goto 800 
      do 303 kk=1,79
      do 304 II=1,57
      J(KK,II)=c1*J(KK,II)
304   continue
303   continue
c      write(6,*) 'c1= ',c1

C calculates atmospheric oxygen, nitrogen and hydrogen values 
c in molecule cm-3 
      O2=0.209460*2.46e19*298./T
      N2 =0.780840*2.46E19*298./T
      H2 =0.5*2.46E13*298./T
      M = O2+N2

c calls subroutines, which calculate deposition velocities (s-1), 
c emission rates (molecule cm-2 s-1) and reaction rate constants
c s-1 for   cm3 molecule-1 s-1 for   cm6 molecule -2 s-1 for
      CALL DEPOSNOP(VD,H,V)
      CALL EMISNOP(EMI1,H,EMI)
      CALL RATECONS(T,H2O,O2,N2,H2,M,K)
      if (par.eq.1) then 
	CALL WASHOUT(T,H2O,WOUT)
      else
       wout=0
      endif
C calls subroutines, which give fotodissosiation rate coefficients (s-1)
      CALL O3(ALFA,J,T,KF(1),KF(2))
      CALL NO2(ALFA,J,KF(3))
      CALL NO3(ALFA,J,KF(4),KF(5))
      CALL N2O5(ALFA,J,T,KF(6))
      CALL H2O2(ALFA,J,KF(7))
      CALL HNO3(ALFA,J,KF(8))
      CALL HCHO(ALFA,J,KF(9),KF(10))
      CALL CH3CHO(ALFA,J,KF(11))
      CALL CH3COC2H5(ALFA,KF(12))
      CALL CH3COCOCH3(ALFA,KF(13))
      CALL CHOCHO(ALFA,J,KF(14))
      CALL CH3OOH(ALFA,J,KF(15))
      CALL CH3COCHO(ALFA,J,KF(16))

      if (i.eq.1) then
	lkm=66
      else 
	lkm=60
      endif
	DO 999 I7=1,lkm
	TOL=1.e-4
	IFAIL = 0
	write(6,'(E10.4)') TT
C calls D02EAF library routine
	CALL D02EAF(TT,TEND,N,Y,TOL,FCN,W,IW,IFAIL)
	DO 233 I4 = 1,54
c        IF (Y(I4).LE.0.) Y(I4)=0.
	IF (Y(I4).LE.1.d-5) Y(I4)=1.d-5
233     CONTINUE
	TT = TEND
	IF (TEND.LT.RAJA) THEN
	TEND = 10.0*TEND
	ELSE
	TEND = TEND+60.
	END IF
      IF (MOD(NINT(TT-3601),3600).EQ.0.and.nint(tt).ne.1) THEN
       WRITE(13,'(28e12.5)') TT-1,(Y(I1),I1=1,27)
       write(14,'(28e12.5)') TT-1,(Y(I1),I1=28,54)
      endif
999   continue
200   CONTINUE
      CLOSE(10)
      CLOSE(11)
      CLOSE(15)
      CLOSE(16)
800   STOP
      END

C***********************************************************************
C Routine for evaluating right hand sides of differential equations.
      SUBROUTINE FCN(TT,Y,F)
      IMPLICIT DOUBLE PRECISION ( A - Z)
C N is the number of the differential equations
      INTEGER N,i3
      PARAMETER (N = 67)
C TT is the time (s)
      DOUBLE PRECISION TT,Y(N),F(N)
      DOUBLE PRECISION K(123),KF(16),H2O,O2,N2,H2,M
      COMMON /LOPUT/ K,KF,H2O,O2,N2,H2,M
      DOUBLE PRECISION V(67),EMI(67),VD(67),EMI1(67),wout
      COMMON /NOP/ V,EMI,VD,EMI1,wout
      DOUBLE PRECISION APU(139)        
      DOUBLE PRECISION T
      COMMON / TEMP /T
      
c      do 256 i3=1,67
c      if (y(i3).lt.1.e-28) y(i3)=1.e-28
c256   continue

C kinetics
      APU(1) = K(1)*Y(3)*O2
      APU(2) = K(2)*Y(3)*Y(5)
      APU(3) = K(3)*Y(2)*O2
      APU(4) = K(4)*H2O*Y(2)
      APU(5) = K(5)*Y(1)*Y(5)
      APU(6) = K(6)*Y(1)*Y(6)
      APU(7) = K(7)*Y(1)*Y(9)
      APU(8) = K(8)*Y(1)*Y(10)
      APU(9) = K(9)*Y(5)*Y(7)
      APU(10) = K(10)*Y(5)*Y(10)
      APU(11) = K(11)*Y(6)*Y(7)
      APU(12) = K(12)*Y(6)*Y(7)
      APU(13) = K(13)*Y(6)*Y(9)
      APU(14) = K(14)*Y(7)*Y(12)
      APU(15) = K(15)*2*Y(7)**2
      APU(16) = K(16)*Y(8)*H2O
      APU(17) = K(17)*Y(8)
      APU(18) = K(18)*Y(9)*Y(10)
      APU(19) = K(19)*Y(9)*Y(12)
      APU(20) = K(20)*Y(9)*H2
      APU(21) = K(21)*Y(9)*Y(13)
      APU(22) = K(22)*2*Y(10)**2
      APU(23) = K(23)*Y(10)*Y(7)
      APU(24) = K(24)*Y(10)*Y(7)
      APU(25) = K(25)*Y(9)*Y(14)
      APU(26) = K(26)*Y(19)*Y(14)
      APU(27) = K(27)*Y(9)*Y(17)
      APU(28) = K(28)*Y(19)*Y(5)
      APU(29) = K(29)*2*Y(19)**2
      APU(30) = K(30)*2*Y(19)**2
      APU(31) = K(31)*Y(10)*Y(19)
      APU(32) = K(32)*Y(9)*Y(20)
      APU(33) = K(33)*Y(7)*Y(20)
      APU(34) = K(34)*Y(9)*Y(15)
      APU(35) = K(35)*Y(21)*Y(9)
      APU(36) = K(36)*Y(22)*Y(5)
      APU(37) = K(37)*Y(23)
      APU(38) = K(38)*Y(23)*O2
      APU(39) = K(39)*Y(22)*Y(19)
      APU(40) = K(40)*Y(9)*Y(24)
      APU(41) = K(41)*Y(25)*Y(6)
      APU(42) = K(42)*Y(26)
      APU(43) = K(43)*Y(25)*Y(5)
      APU(44) = K(44)*Y(19)*Y(25)
      APU(45) = K(45)*2*Y(22)**2
      APU(46) = K(46)*2*Y(25)**2
      APU(47) = K(47)*Y(27)*Y(9)
      APU(48) = K(48)*Y(28)*Y(5)
      APU(49) = K(49)*Y(29)*O2
      APU(50) = K(50)*Y(29)
      APU(51) = K(51)*Y(28)*Y(19)
      APU(52) = K(52)*Y(30)*Y(9)
      APU(53) = K(53)*Y(31)*Y(5)
      APU(54) = K(54)*Y(31)*Y(19)
      APU(55) = K(55)*Y(34)*Y(9)
      APU(56) = K(56)*Y(33)*Y(5)
      APU(57) = K(57)*Y(33)*Y(19)
      APU(58) = K(58)*Y(34)*Y(1)
      APU(59) = K(59)*Y(35)*Y(1)
      APU(60) = K(60)*Y(35)*Y(1)
      APU(61) = K(61)*Y(35)*Y(9)
      APU(62) = K(62)*Y(36)*Y(5)
      APU(63) = K(63)*Y(36)*Y(19)
      APU(64) = K(64)*Y(9)*Y(37)
      APU(65) = K(65)*Y(47)*Y(5)
      APU(66) = K(66)*Y(38)*Y(9)
      APU(67) = K(67)*Y(39)*Y(5)
      APU(68) = K(68)*Y(9)*Y(41)
      APU(69) = K(69)*Y(9)*Y(40)
      APU(70) = K(70)*Y(9)*Y(42)
      APU(71) = K(71)*Y(43)*Y(5)
      APU(72) = K(72)*Y(9)*Y(44)
      APU(73) = K(73)*Y(45)*Y(5)
      APU(74) = KF(1)*Y(1)
      APU(75) = KF(2)*Y(1)
      APU(76) = KF(3)*Y(6)
      APU(77) = KF(4)*Y(7)
      APU(78) = KF(5)*Y(7)
      APU(79) = KF(6)*Y(8)
      APU(80) = KF(7)*Y(12)
      APU(81) = KF(8)*Y(13)
      APU(82) = KF(9)*Y(20)
      APU(83) = KF(10)*Y(20)
      APU(84) = KF(11)*Y(24)
      APU(85) = KF(12)*Y(30)
      APU(86) = KF(13)*Y(32)
      APU(87) = KF(14)*Y(41)
      APU(88) = KF(15)*Y(46)
      APU(89) = KF(16)*Y(40)
      APU(90) = K(74)*Y(18)*O2
      APU(91) = K(75)*Y(48)*O2
      APU(92) = K(76)*Y(50)*O2
      APU(93) = K(77)*Y(49)*H2O
      APU(94) = K(78)*Y(42)*Y(1)
      APU(95) = K(79)*Y(42)*Y(7)
      APU(96) = K(80)*Y(51)*Y(9)
      APU(97) = K(81)*Y(51)*Y(7)
      APU(98) = K(82)*Y(51)*Y(1)
      APU(99) = K(83)*Y(52)*Y(9)
      APU(100) = K(84)*Y(52)*Y(7)
      APU(101) = K(85)*Y(52)*Y(1)
      APU(102) = K(86)*Y(53)*Y(9)
      APU(103) = K(87)*Y(53)*Y(7)
      APU(104) = K(88)*Y(53)*Y(1)
      APU(105) = K(89)*Y(54)*Y(9)
      APU(106) = K(90)*Y(54)*Y(7)
      APU(107) = K(91)*Y(54)*Y(1)
      APU(108) = K(92)*Y(55)*Y(9)
      APU(109) = K(93)*Y(55)*Y(9)
      APU(110) = K(94)*Y(57)
      APU(111) = K(95)*Y(56)*Y(6)
      APU(112) = K(96)*Y(56)*Y(1)
      APU(113) = K(97)*Y(56)*O2
      APU(114) = K(98)*Y(59)
      APU(115) = K(99)*Y(61)*Y(6)
      APU(116) = K(100)*Y(64)
      APU(117) = K(101)*Y(59)
      APU(118) = K(102)*Y(59)*Y(5)
      APU(119) = K(103)*Y(58)*Y(1)
      APU(120) = K(104)*Y(58)*Y(6)
      APU(121) = K(105)*Y(58)*Y(6)
      APU(122) = K(106)*Y(58)*O2
      APU(123) = K(107)*Y(61)
      APU(124) = K(108)*Y(61)*Y(5)
      APU(125) = K(109)*Y(60)*Y(6)
      APU(126) = K(110)*Y(60)*Y(1)
      APU(127) = K(111)*Y(60)*O2
      APU(128) = K(112)*Y(63)
      APU(129) = K(113)*Y(63)*Y(6)
      APU(130) = K(114)*Y(65)
      APU(131) = K(115)*Y(60)
      APU(132) = K(116)*Y(62)
      APU(133) = K(117)*Y(62)
      APU(134) = K(118)*Y(63)*Y(5)
      APU(135) = K(119)*Y(59)*Y(1)
      APU(136) = K(120)*Y(59)*Y(6)
      APU(137) = K(121)*Y(66)
      APU(138) = K(122)*Y(56)*O2
      APU(139) = K(123)*Y(59)*O2

C Differential equations for concentrations
c O3
      F(1) = APU(1)-APU(5)-APU(6)-APU(7)-APU(8)-APU(58)-APU(59)
     1       -APU(60)-APU(74)-APU(75)-APU(94)-APU(98)-APU(101)
     2       -APU(104)-APU(107)-APU(112)-APU(119)-APU(126)
     2       -APU(135)-V(1)*Y(1)
c O(D)
      F(2) = -APU(3)-APU(4)+APU(74)-V(2)*Y(2)
c O(P)
      F(3) = -APU(1)-APU(2)+APU(3)+APU(75)+APU(76)+APU(77)-V(3)*Y(3)
c O2
      F(4)=0
c NO
      F(5) = -APU(2)-APU(5)-APU(9)-APU(10)+APU(11)-APU(28)-APU(36)
     1       -APU(43)-APU(48)-APU(53)-APU(56)-APU(62)-APU(65)-APU(67)
     2       -APU(71)-APU(73)+APU(76)+APU(78)-apu(96)-apu(99)-apu(102)
     3       -apu(105)+APU(111)-APU(118)+APU(120)+APU(121)-APU(124)
     4       +APU(125)-APU(134)-V(5)*Y(5)+emi(5)
c NO2
      F(6) = APU(2)+APU(5)-APU(6)-APU(12)-APU(13)+2*APU(9)+APU(10)
     1       +2*APU(15)+APU(17)+APU(24)+APU(28)+APU(36)-APU(41)+APU(42)
     2       +APU(43)+APU(48)+APU(53)+APU(56)+APU(62)+APU(65)+APU(67)
     3       +APU(71)+APU(73)-APU(76)+APU(77)+APU(79)+APU(81)+apu(96)
     4       +apu(99)+apu(102)+apu(105)-APU(111)-APU(115)+APU(116)
     5       +APU(118)-APU(120)-APU(121)+APU(124)-APU(125)-APU(129)
     6       +APU(130)+APU(134)-APU(136)+APU(137)-V(6)*Y(6)+EMI(6)
c NO3
      F(7) = APU(6)-APU(9)-APU(11)-APU(12)-APU(14)-2*APU(15)+APU(17)
     1       +APU(21)-APU(23)-APU(24)-APU(33)-APU(77)-APU(78)+APU(79)
     2       -APU(95)-APU(97)-APU(100)-APU(103)-APU(106)-V(7)*Y(7)
c N2O5
      F(8) = APU(12)-APU(16)-APU(17)-APU(79)-V(8)*Y(8)
c OH
      F(9) = 2*APU(4)-APU(7)+APU(8)+APU(10)-APU(13)-APU(18)-APU(19)
     1       -APU(20)-APU(21)+APU(24)-APU(25)-APU(27)-APU(32)-APU(34)
     2       -APU(35)-APU(40)-APU(47)-APU(52)-APU(55)+0.19*APU(59)
     3       -APU(61)-APU(64)-apu(66)-APU(68)-APU(69)-APU(70)-APU(72)
     4       +2*APU(80)+APU(81)+APU(88)-APU(108)-APU(109)-APU(96)
     5       -APU(102)-APU(105)-apu(99)-V(9)*Y(9)
c HO2
      F(10) = APU(7)-APU(8)-APU(10)+APU(14)-APU(18)+APU(19)-2*APU(22)
     1       -APU(23)-APU(24)+APU(91)-APU(31)+apu(32)+APU(33)
     2       +apu(34)+APU(38)+APU(49)+APU(51)+APU(53)+2*APU(54)
     3       +APU(56)+0.12*APU(58)+0.29*APU(59)+0.12*APU(60)
     4       +APU(62)+APU(63)+APU(65)+APU(67)+APU(71)+APU(73)
     5       +2*APU(82)+APU(84)+APU(89)+APU(92)+apu(96)+apu(99)
     6       +apu(102)+apu(105)-V(10)*Y(10)
c H2SO4
      F(11) = APU(93)-V(11)*Y(11)-wout*Y(11)
c H2O2
      F(12) = -APU(14)-APU(19)+APU(22)-APU(80)-V(12)*Y(12)-wout*Y(12)
c HNO3
      F(13) = APU(13)+APU(14)+2*APU(16)-APU(21)+APU(23)+APU(33)
     1        -APU(81)-V(13)*Y(13)-wout*Y(13)
c SO2
      F(14) = -APU(25)-APU(26)+APU(121)+APU(131)+APU(138)+APU(139)
     1         -V(14)*Y(14)+EMI(14)
c       F(14) = 0
c CO
      F(15) = APU(32)+APU(33)-APU(34)+0.42*APU(58)+0.24*APU(59)
     1        +0.42*APU(60)+2*APU(68)+APU(69)+APU(82)+APU(83)+APU(84)
     2        +APU(87)+APU(89)-V(15)*Y(15)+emi(15)
c CO2
      F(16)=0
c CH4
      F(17)=0
c CH3
      F(18) = APU(27)-APU(90)+APU(37)+APU(43)+APU(44)+2*APU(46)
     1        +APU(138)+APU(139)+APU(131)+APU(132)-V(18)*Y(18)
c CH3O2
      F(19) = -APU(26)+APU(90)-APU(28)-2*APU(29)-2*APU(30)-APU(31)
     1        -APU(39)-APU(44)-APU(51)-APU(54)+0.43*APU(59)-APU(63)
     2        +APU(84)-APU(57)-V(19)*Y(19)
c HCHO
      F(20) = APU(91)+APU(30)-APU(32)-APU(33)+APU(37)+APU(51)+APU(54)
     1        +2*APU(56)+APU(58)+APU(59)+APU(62)+APU(63)
     2        +APU(71)+APU(73)-APU(82)-APU(83)+APU(87)
     3        -V(20)*Y(20)+EMI(20)
c C2H6
      F(21) = -APU(35)+EMI(21)-V(21)*Y(21)
c C2H5O2
      F(22) = APU(35)-APU(36)-APU(39)-2*APU(45)+APU(85)-V(22)*Y(22)
c C2H5O
      F(23) = APU(36)-APU(37)-APU(38)+APU(39)+2*APU(45)-V(23)*Y(23)
c CH3CHO
      F(24) = APU(38)-APU(40)+APU(50)+APU(60)+APU(62)-APU(84)+APU(63)
     1        -V(24)*Y(24)+emi(24)
c CH3COO2
      F(25) = APU(40)-APU(41)+APU(42)-APU(43)-APU(44)-2*APU(46)+APU(69)
     1        +APU(85)+2*APU(86)-V(25)*Y(25)
c PAN
      F(26) = APU(41)-APU(42)-V(26)*Y(26)
c nC4H10
      F(27) = -APU(47)+EMI(27)-V(27)*Y(27)
c C4H9O2
      F(28) = APU(47)-APU(48)-APU(51)-V(28)*Y(28)
c C4H9O
      F(29) = APU(48)-APU(49)-APU(50)+APU(51)-V(29)*Y(29)
c CH3COC2H5
      F(30) = APU(49)-APU(52)-APU(85)-V(30)*Y(30)+emi(30)
c CH3COCHO2CH3
      F(31) = APU(52)-APU(53)-APU(54)-V(31)*Y(31)
c CH3COCOCH3
      F(32) = APU(53)+APU(54)-APU(86)-V(32)*Y(32)
c CH2O2CH2OH
      F(33) = APU(55)-APU(56)-APU(57)-V(33)*Y(33)
c C2H4
      F(34) = -APU(55)-APU(58)+EMI(34)-v(34)*Y(34)
c C3H6
      F(35) = -APU(59)-APU(60)-APU(61)+EMI(35)-V(35)*Y(35)
c CH3CHO2CH2OH
      F(36) = APU(61)-APU(62)-APU(63)-V(36)*Y(36)
c o-ksyleeni
      F(37) = -APU(64)+EMI(37)-V(37)*Y(37)
c CHOCH=CHCOCH3
      F(38) = APU(65)-APU(66)-V(38)*Y(38)
c CH3COCHOH-CHO2CHO
      F(39) = APU(66)-APU(67)-V(39)*Y(39)
c CH3COCHO
      F(40) = APU(73)+APU(65)-APU(69)-APU(89)+APU(67)-V(40)*Y(40)
c CHOCHO
      F(41) = APU(67)-APU(68)-APU(87)-V(41)*Y(41)
c C5H8
      F(42) = -APU(70)-APU(94)-APU(95)+EMI(42)-V(42)*Y(42)
c OHC5H8O2
      F(43) = APU(70)-APU(71)-V(43)*Y(43)
c CH3COCH=CH2
      F(44) = APU(71)-APU(72)-V(44)*Y(44)
c OHCH3COCHCH2O2
      F(45) = APU(72)-APU(73)-V(45)*Y(45)
c CH3OOH
      F(46) = APU(31)-APU(88)-V(46)*Y(46)-wout*Y(46)
c
      F(47) = APU(64)-APU(65)
c CH3O
      F(48) = APU(26)+APU(28)+0.05*APU(59)+2*APU(29)+APU(39)+APU(44)
     1        +APU(57)+APU(63)+apu(88)-APU(91)-V(48)*Y(48)
C SO3
      F(49) = APU(26)+APU(92)-APU(93)+APU(132)-V(49)*Y(49)
C HSO3
      F(50) = APU(25)-APU(92)-V(50)*Y(50)
C ALFAPINENE
      F(51) = -APU(96)-APU(97)-APU(98)+EMI(51)-V(51)*Y(51)
C BEETAPINENE
      F(52) = -APU(99)-APU(100)-APU(101)+EMI(52)-V(52)*Y(52)
C 3-KARENE
      F(53) = -APU(102)-APU(103)-APU(104)+EMI(53)-V(53)*Y(53)
C D-LIMONENE
      F(54) = -APU(105)-APU(106)-APU(107)+EMI(54)-V(54)*Y(54)
C DMS
c      F(55) = -APU(108)-APU(109)+EMI(55)-V(55)*Y(55)
       f(55)=0
C CH3S
c      F(56) = APU(108)-APU(111)-APU(112)-APU(113)+APU(114)-APU(138)
c     1       -V(56)*Y(56)
      f(56)=0
C CH3S(OH)CH3
c      F(57) = APU(109)-APU(110)-V(57)*Y(57)
      f(57)=0
C CH3SO
c      F(58) = APU(110)+APU(111)+APU(112)+APU(118)-APU(119)-APU(120)
c     1        -APU(121)-APU(122)+APU(123)+APU(135)-V(58)*Y(58)
      f(58)=0
C CH3SOO
c      F(59) = APU(113)-APU(114)-APU(117)-APU(118)-APU(135)-APU(136)
c     1        +APU(137)-APU(139)-V(59)*Y(59)
      f(59)=0
C CH3SO2
c      F(60) = APU(117)+APU(119)+APU(120)+APU(124)-APU(125)-APU(126)
c     1        -APU(127)+APU(128)-V(60)*Y(60)
      f(60)=0
C CH3S(O)O2
c     F(61) = -APU(115)+APU(116)+APU(122)-APU(123)-APU(124)-V(61)*Y(61)
      f(61)=0
C CH3SO3
c      F(62) = APU(125)+APU(126)-APU(132)-APU(133)+APU(134)-V(62)*Y(62)
      f(62)=0
C CH3S(O)2O2
c      F(63) = APU(127)-APU(128)-APU(129)+APU(130)-APU(134)-V(63)*Y(63)
      f(63)=0
C CH3S(O)O2NO2
c      F(64) = APU(115)-APU(116)-V(64)*Y(64)
      f(64)=0
C CH2S(O)2O2NO2
c      F(65) = APU(129)-APU(130)-V(65)*Y(65)
      f(65)=0
C CH3SOONO2
c      F(66) = APU(136)-APU(137)-V(66)*Y(66)
      f(66)=0
C CH3SO3H (MSA)
c      F(67) = APU(133)-V(67)*Y(67)
      f(67)=0
c      do 257 i3=1,54      
c      if (f(i3).gt.1.e20) write(6,*) f(i3)
c257   continue
      RETURN
      END
C

C subroutines for the rate constants of the photodissociation reactions
      SUBROUTINE O3(ALFA,J,T,K1,K2)
      IMPLICIT DOUBLE PRECISION (A - Z)
      INTEGER I,L
      DOUBLE PRECISION J(79,57),ALFA
      DOUBLE PRECISION T
      K1=0.
      DL=5.
      VA=1.E-20*1.E14
      T1=T-230.
      A=0.332+2.565E-4*T1+1.152E-5*T1**2+2.313E-8*T1**3
      B=-0.575+5.59E-3*T1-1.439E-5*T1**2+3.27E-8*T1**3
      C=0.466+8.883E-4*T1-3.546E-5*T1**2-3.519E-7*T1**3
      LO=308.20+4.4871E-2*T1+6.9380E-5*T1*2-2.5452E-6*T1**3
      DO 500 I=1,7
      READ(20,*) AL,CR
      L=NINT(AL)
      IF (L.LE.300) THEN
	FI=0.9
      ELSE
	FI=A*ATAN(B*(REAL(L)-LO))+C
      END IF
      IF (FI.GT.0.9) FI=0.9
      IF (FI.LT.0.) FI=0.
      IF (ALFA.GT.86.) THEN
       K1=K1+CR*FI*J((L-290)/5+1,57)*DL*VA
       K1 = -K1/4*(ALFA-90.)
      else
       K1=K1+CR*FI*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA 
      endif
500   CONTINUE
      CLOSE(20)
      IF (ALFA.GT.89.9) THEN
       KULMA=89.9
       K2=1.23E-3*DEXP(-0.6/COS(KULMA*3.14159265/180.))
      ELSE 
       K2=1.23E-3*DEXP(-0.6/COS(ALFA*3.14159265/180.))
      ENDIF
      if (alfa.ge.90.) k2=0.
      RETURN
      END
C
      SUBROUTINE NO2(ALFA,J,K)
      IMPLICIT DOUBLE PRECISION (A - Z)
      INTEGER I,L
      DOUBLE PRECISION J(79,57),ALFA
      K=0.
      DL=5.
      VA=1.E-20*1.E14
      DO 500 I=1,25
      READ(21,*) L,CR,FI
      IF (ALFA.GT.86.) THEN
       K=K+CR*FI*J((L-290)/5+1,57)*DL*VA
       K = -K/4*(ALFA-90.)
      ELSE
       K=K+CR*FI*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
      ENDIF
  500 CONTINUE
      CLOSE(21)
      RETURN
      END
C
      SUBROUTINE NO3(ALFA,J,K1,K2)
      IMPLICIT DOUBLE PRECISION (A - Z)
      INTEGER I,L
      DOUBLE PRECISION J(79,57),K1,K2
      K1=0.
      K2=0.
      DL=5.
      VA=1.E-19*1.E14
      DO 500 I=1,32
      READ(22,*) AL,TULO1,TULO2
      L=NINT(AL)
      IF (ALFA.GT.86.) THEN
       K1=K1+TULO1*J((L-290)/5+1,57)*DL*VA
       K1 = -K1/4*(ALFA-90.)
       K2=K2+TULO2*J((L-290)/5+1,57)*DL*VA
       K2 = -K2/4*(ALFA-90.)
      ELSE
       K1=K1+TULO1*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
       K2=K2+TULO2*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
      ENDIF
  500 CONTINUE
      CLOSE(22)
      RETURN
      END
C
      SUBROUTINE N2O5(ALFA,J,T,K)
      IMPLICIT DOUBLE PRECISION (A - Z)
      INTEGER I,L
      DOUBLE PRECISION J(79,57),ALFA
      K=0.
      DL=5.
      L=290
      FI=1.
      VA=1.E-20*1.E14
      DO 500 I=1,13
      CR=DEXP(2.735+(4728-17.13*L)/T)
      IF (ALFA.GT.86.) THEN
       K=K+CR*FI*J((L-290)/5+1,57)*DL*VA
       K = -K/4*(ALFA-90.)
      ELSE
       K=K+CR*FI*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
      ENDIF
      L=L+5
  500 CONTINUE
      RETURN
      END
C
      SUBROUTINE HNO3(ALFA,J,K)
      IMPLICIT DOUBLE PRECISION (A - Z)
      INTEGER I,L
      DOUBLE PRECISION J(79,57),ALFA
      K=0.
      DL=5.
      VA=1.E-20*1.E14
      DO 500 I=1,7
      READ(24,*) L,CR,FI
      IF (ALFA.GT.86.) THEN
       K=K+CR*FI*J((L-290)/5+1,57)*DL*VA
       K = -K/4*(ALFA-90.)
      ELSE
       K=K+CR*FI*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
      ENDIF
  500 CONTINUE
      CLOSE(24)
      RETURN
      END
C
      SUBROUTINE H2O2(ALFA,J,K)
      IMPLICIT DOUBLE PRECISION (A - Z)
      INTEGER I,L
      DOUBLE PRECISION J(79,57),ALFA
      K=0.
      DL=5.
      VA=1.E-20*1.E14
      DO 500 I=1,13
      READ(25,*) L,CR,FI
      IF (ALFA.GT.86.) THEN
       K=K+CR*FI*J((L-290)/5+1,57)*DL*VA
       K = -K/4*(ALFA-90.)
      ELSE
       K=K+CR*FI*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
      ENDIF
  500 CONTINUE
      CLOSE(25)
      RETURN
      END
C
      SUBROUTINE HCHO(ALFA,J,K1,K2)
      IMPLICIT DOUBLE PRECISION (A - Z)
      INTEGER I,L
      DOUBLE PRECISION J(79,57),ALFA
      K1=0.
      K2=0.
      DL=5.
      VA=1.E-20*1.E14
      DO 500 I=1,15
      READ(26,*) AL,CR,FI1,FI2
      L=NINT(AL)
      IF (ALFA.GT.86.) THEN
       K1=K1+CR*FI1*J((L-290)/5+1,57)*DL*VA
       K2=K2+CR*FI2*J((L-290)/5+1,57)*DL*VA
       K1 = -K1/4*(ALFA-90.)
       K2 = -K2/4*(ALFA-90.)
      ELSE
       K1=K1+CR*FI1*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
       K2=K2+CR*FI2*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
      ENDIF
  500 CONTINUE
      CLOSE(26)
      RETURN
      END
C
      SUBROUTINE CH3CHO(ALFA,J,K)
      IMPLICIT DOUBLE PRECISION (A - Z)
      INTEGER I,L
      DOUBLE PRECISION J(79,57),ALFA
      K=0.
      DL=5.
      VA=1.E-20*1.E14
      DO 500 I=1,8
      READ(27,*) L,CR,FI
      IF (ALFA.GT.86.) THEN
       K=K+CR*FI*J((L-290)/5+1,57)*DL*VA
       K = -K/4*(ALFA-90.)
      ELSE
       K=K+CR*FI*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
      ENDIF
  500 CONTINUE
      CLOSE(27)
      RETURN
      END
C
      SUBROUTINE CHOCHO(ALFA,J,K)
      IMPLICIT DOUBLE PRECISION (A - Z)
      INTEGER I,L
      DOUBLE PRECISION J(79,57),ALFA
      K=0.
      DL=5.
      VA=1.E-20*1.E14
      DO 500 I=1,35
      READ(28,*) AL,CR
      L=NINT(AL)
      IF (L.LT.325) THEN
	FI=0.4
      ELSE
	FI=0.029
      ENDIF
      IF (ALFA.GT.86.) THEN
	K=K+CR*FI*J((L-290)/5+1,57)*DL*VA
	K = -K/4*(ALFA-90.)
      ELSE
	K=K+CR*FI*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
      ENDIF
500   CONTINUE
      CLOSE(28)
      RETURN
      END
C
      SUBROUTINE CH3OOH(ALFA,J,K)
      IMPLICIT DOUBLE PRECISION (A - Z)
      INTEGER I,L
      DOUBLE PRECISION J(79,57),ALFA
      K=0.
      DL=5.
      VA=1.E-20*1.E14
      DO 500 I=1,13
      READ(29,*) AL,CR,FI
      L=NINT(AL)
      IF (ALFA.GT.86.) THEN
       K=K+CR*FI*J((L-290)/5+1,57)*DL*VA
       K = -K/4*(ALFA-90.)
      ELSE
	K=K+CR*FI*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
      ENDIF
  500 CONTINUE
      CLOSE(29)
      RETURN
      END
C
      SUBROUTINE CH3COCHO(ALFA,J,K)
      IMPLICIT DOUBLE PRECISION (A - Z)
      INTEGER I,L
      DOUBLE PRECISION J(79,57),ALFA
      K=0.
      DL=5.
      VA=1.E-20*1.E14
      FI=0.107
      DO 500 I=1,29
      READ(30,*) AL,CR
      L=NINT(AL)
      IF (ALFA.GT.86.) THEN
       K=K+CR*FI*J((L-290)/5+1,57)*DL*VA
       K = -K/4*(ALFA-90.)
      ELSE
       K=K+CR*FI*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
      ENDIF
  500 CONTINUE
      CLOSE(30)
      RETURN
      END
C
      SUBROUTINE CH3COC2H5(ALFA,K)
      IMPLICIT DOUBLE PRECISION (A - Z)
      DOUBLE PRECISION ALFA
      A=2.43E-5
      B=0.877
      IF (ALFA.GT.89.9) THEN
       K=A*DEXP(-B/COS(89.9*3.14159265/180.))
      ELSE
       K=A*DEXP(-B/COS(ALFA*3.14159265/180.))
      ENDIF
      if (alfa.ge.90.) k=0.
      RETURN
      END
C
      SUBROUTINE CH3COCOCH3(ALFA,K)
      IMPLICIT DOUBLE PRECISION (A - Z)
      DOUBLE PRECISION ALFA
      A=2.7E-4
      B=0.79
      IF (ALFA.GT.89.9) THEN
       K=A*DEXP(-B/COS(89.9*3.14159265/180.))
      ELSE
       K=A*DEXP(-B/COS(ALFA*3.14159265/180.))
      ENDIF
      if (alfa.ge.90.) k=0.
      RETURN
      END

C Rate constants
      SUBROUTINE RATECONS(T,H2O,O2,N2,H2,M,K)
      IMPLICIT DOUBLE PRECISION (A - Z)
      DOUBLE PRECISION T,K(123),KF(16),O2,N2,M

	KoO2=6.2E-34*(T/300.)**(-2.8)*O2
	KoN2=5.6E-34*(T/300.)**(-2.8)*N2
	KINF=2.8E-12
	Fc=DEXP(-T/696.)
	KO2=(KoO2/(1.+KoO2/KINF))*Fc**(1/(1+(DLOG10(KoO2/KINF))**2))
	KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
      K(1) = KO2+KN2
	KoO2=8.6E-32*(T/300.)**(-1.8)*O2
	KoN2=1.0E-31*(T/300.)**(-1.6)*N2
	KINF=3.0E-11*(T/300.)**0.3
	Fc=DEXP(-T/1850.)
	KO2=(KoO2/(1.+KoO2/KINF))*Fc**(1/(1+(DLOG10(KoO2/KINF))**2))
	KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
      K(2) = KO2+KN2
      K(3) = 3.2E-11*DEXP(67./T)
      K(4) = 2.2E-10
      K(5) = 1.8E-12*DEXP(-1370./T)
      K(6) = 1.2E-13*DEXP(-2450./T)
      K(7) = 1.9E-12*DEXP(-1000./T)
      K(8) = 1.4E-14*DEXP(-600./T)
      K(9) = 1.8E-11*DEXP(110./T)
      K(10) = 3.7E-12*DEXP(240./T)
      K(11) = 2.3E-12*DEXP(-1000./T)
	KoN2=2.7E-30*(T/300.)**(-3.4)*N2
	KINF=2.0E-12*(T/300.)**0.2
	Fc = DEXP(-T/250.)+DEXP(-1050./T)
	KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
      K(12) = KN2
	KoO2=2.2E-30*(T/300.)**(-2.9)*O2
	KoN2=2.6E-30*(T/300.)**(-2.9)*N2
	KINF=5.2E-11
	Fc=DEXP(-T/353.)
	KO2=(KoO2/(1.+KoO2/KINF))*Fc**(1/(1+(DLOG10(KoO2/KINF))**2))
	KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
      K(13) = KO2+KN2
      K(13) = KN2
      K(14) = 4.1E-16
      K(15) = 8.5E-13*DEXP(-2450./T)
      K(16) = 1.3E-21
	KoN2=2.2E-3*(T/300.)**(-4.4)*DEXP(-11080./T)*N2
	KINF=9.7E14*(T/300.)**0.1*DEXP(-11080./T)
	Fc=DEXP(-T/250.)+DEXP(-1050./T)
	KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
      K(17) = KN2
      K(18) = 4.8E-11*DEXP(250./T)
      K(19) = 2.9E-12*DEXP(-160./T)
      K(20) = 7.7E-12*DEXP(-2100./T)
      K(21) = 9.4E-15*DEXP(778./T)
      K(22) = 1.6E-12
      K(23) = 4.3E-12
      K(24) = 4.3E-12
	KoO2=5.0E-31*(T/300.)**(-3.3)*O2
	KoN2=5.0E-31*(T/300.)**(-3.3)*N2
	KINF=2.E-12
	Fc=DEXP(-T/380)
	KO2=(KoO2/(1.+KoO2/KINF))*Fc**(1/(1+(DLOG10(KoO2/KINF))**2))
	KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
      K(25) = KO2+KN2
      K(25) = KO2
      K(26) = 4.0E-17
      K(27) = 3.9E-12*DEXP(-1885./T)
      K(28) = 4.2E-12*DEXP(180./T)
      K(29) = 1.1E-13*DEXP(365./T)
      K(30) = 1.1E-13*DEXP(365./T)
      K(31) = 3.8E-13*DEXP(780./T)
      K(32) = 8.8E-12*DEXP(25./T)
      K(33) = 5.8E-16
      K(34) = 2.4E-13
      K(35) = 7.8E-12*DEXP(-1020./T)
      K(36) = 8.9e-12
      K(37) = 33.0
      K(38) = 9.5E-15
      K(39) = 2.5E-14
      K(40) = 5.6E-12*DEXP(310./T)
	KoN2=2.7E-28*(T/300.)**(-7.1)*N2
	KINF=1.2E-11*(T/300.)**(-0.9)
	Fc=0.3
	KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
      K(41) = KN2
	KoN2=4.9E-3*DEXP(-12100./T)*N2
	KoN2=4.9E-3*DEXP(-12100./T)*M
	KINF=4.0E16*DEXP(-13600./T)
	Fc=0.3
	KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
      K(42) = KN2
      K(43) = 2.0E-11
      K(44) = 5.5E-12
      K(45) = 9.8E-14*DEXP(-110./T)
      K(46) = 2.8E-12*DEXP(530./T)
      K(47) = 1.4E-11*DEXP(-559./T)
      K(48) = 3.0E-12
      K(49) = 2.1E-16
      K(50) = 1.2E3
      K(51) = 2.5E-14
      K(52) = 8.8E-13
      K(53) = 3.1E-13
      K(54) = 2.5E-14
	KoO2=9.5E-29*(T/300.)**(-3.1)*O2
	KoN2=9.5E-29*(T/300.)**(-3.1)*N2
	KINF=9.E-12
	Fc=DEXP(-T/840.)
	KO2=(KoO2/(1.+KoO2/KINF))*Fc**(1/(1+(DLOG10(KoO2/KINF))**2))
	KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
      K(55) = KO2+KN2
      K(55) = KN2
      K(56) = 3.1E-13
      K(57) = 2.5E-14
      K(58) = 1.2E-14*DEXP(-2630./T)
      K(59) = 6.5E-15*DEXP(-1880./T)
      K(60) = 6.5E-15*DEXP(-1880./T)
	KoO2=8.E-27*(T/300.)**(-3.5)*O2
	KoN2=8.E-27*(T/300.)**(-3.5)*N2
	KINF=3.0E-11
	Fc=DEXP(-T/433.)
	KO2=(KoO2/(1.+KoO2/KINF))*Fc**(1/(1+(DLOG10(KoO2/KINF))**2))
	KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
      K(61) = KO2+KN2
      K(62) = 3.1E-13
      K(63) = 2.5E-14
      K(64) = 1.1E-11
      K(65) = 3.1E-13
      K(66) = 2.0E-11
      K(67) = 3.1E-13
      K(68) = 1.1E-11
      K(69) = 1.7E-11
      K(70) = 2.54E-11*DEXP(410./T)
      K(71) = 3.0E-13
      K(72) = 2.0E-11
      K(73) = 3.0E-13
	KoN2=1.0E-30*(T/300.)**(-3.3)*N2
	KINF=2.2E-12*(T/300.)**1.0
	Fc=0.27
	KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
      K(74) = KN2
      K(75) = 1.9E-15
      K(76) = 4.4E-13
      K(77) = 9.1E-13
      K(78) = 1.2e-17
      K(79) = 3.2e-13
      K(80) = 0.98e-11*dexp(500./T)
      K(81) = 5.8e-12
      K(82) = 4.6e-15*dexp(-1170./T)
      K(83) = 7.95e-11
      K(84) = 2.36e-12
      K(85) = 2.1e-17
      K(86) = 1.6e-11*dexp(500./T)
      K(87) = 10.1e-12
      K(88) = 6.5e-15*dexp(-1170./T)
      K(89) = 16.9e-11
      K(90) = 1.4e-11
      K(91) = 6.4e-16
      K(92) = 1.1D-11*dexp(-240./T)
      K(93) = 1.2D-12
      K(94) = 0.5
      K(95) = 5.6D-11
      K(96) = 5.4D-12
      K(97) = 6.1D-19
      K(98) = 1.
      K(99) = 10**(-11.23)
      K(100) = 10.**(-1.951)     
      K(101) = 5.0
      K(102) = 1.11D-11
      K(103) = 3.D-13
      K(104) = 1.2D-11
      K(105) = 0.85D-11
      K(106) = 7.7D-18
      K(107) = 1.7D2
      K(108) = 10.**(-10.62)
      K(109) = 10.**(-10.92)
      K(110) = 10.**(-12.22)
      K(111) = 2.6D-18
      K(112) = 3.3
      K(113) = 10.**(-11.23)
      K(114) = 10.**(-1.951)
      K(115) = 5.00135
      K(116) = 5.00135
      K(117) = 5.005
      K(118) = 10.**(-10.62)
      K(119) = 10.**(-12.1)
      K(120) = 10.**(-11.30)
      K(121) = 10.**(0.2553)
      K(122) = 3.0D-18
      K(123) = 3.0D-18
      RETURN
      END


      SUBROUTINE DEPOSNOP(VD,H,V)
      IMPLICIT DOUBLE PRECISION (A - Z)
      INTEGER I
      DOUBLE PRECISION VD(67),V(67)
      DOUBLE PRECISION J(79,57),ALFA
      COMMON /INTENS/ J,ALFA

      IF (ALFA.GE.86.) THEN
       V(1)=0.5*VD(1)
      ELSE
       V(1)=VD(1)
      ENDIF
      DO 314 I=2,67
      IF (ALFA.GE.86.) THEN
       V(I)=0.25*VD(I)
      ELSE
       V(I)=VD(I)
      ENDIF
314   CONTINUE

C Deposition velocities divided by mixing height H => s-1
      DO 315 I=1,67
      V(I)=V(I)/H
315   CONTINUE
      RETURN
      END

      SUBROUTINE EMISNOP(EMI1,H,EMI)
      IMPLICIT DOUBLE PRECISION (A - Z)
      INTEGER I
      DOUBLE PRECISION EMI1(67),EMI(67)
      DOUBLE PRECISION J(79,57),ALFA
      COMMON /INTENS/ J,ALFA

      DO 316 I=1,67
       IF (ALFA.LT.86.) THEN
	EMI(I)=1.5*EMI1(I)
       ELSE
	EMI(I)=0.5*EMI1(I)
	EMI(42)=0
       ENDIF
316   CONTINUE

C Emission rates divided by mixing height H => molecule cm-3 s-1
      DO 317 I=1,67
	EMI(I)=EMI(I)/H
317   CONTINUE
      RETURN
      END

      SUBROUTINE washout(T,H2O,WOUT)
c calculates the washout coefficient (s-1)
      implicit double precision (a - z)
      psat = DEXP(77.344913-7235.4247/T-8.2*DLOG(T)+0.0057113*T)
      rh=100.*H2O*T/psat*1.381d-17
      if (rh.gt.100.) rh=100.
      if (rh.ge.95.) then
	wout=1.d-4
      else 
	wout=0.
      endif
c      write(6,*) 'real RH= ',rh, 'wout= ',wout
      return
      END

      SUBROUTINE CLOUDINESS(CL,NHcl,ALFA,C1,koodi)
      IMPLICIT DOUBLE PRECISION (A-Z)
      integer cl,nhcl,koodi
      koodi=0
      if (alfa.gt.86.) then 
       alfa1=86.
      else
       alfa1=alfa
      endif
      M1=950./(1013.*COS(alfa1*3.14159265/180.))
c      write(6,*) 'm1= ',m1
      IF (CL.EQ.0) THEN
	   C1=1.
      ELSE IF (CL.EQ.1) THEN
	   T1=0.2684-0.0101*M1
	   C1=1.-NHCL/8.*(1.-T1)
      ELSE IF (CL.EQ.2.OR.CL.EQ.3) THEN
	   T1=0.3658-0.0149*M1
	   C1=1.-NHCL/8.*(1.-T1)
      ELSE IF (CL.EQ.4) THEN
	   T1=0.2363+0.0145*M1
	   C1=1.-NHCL/8.*(1.-T1)
      ELSE IF (CL.EQ.5) THEN
	   T1=0.4130-0.0014*M1
	   C1=1.-NHCL/8.*(1.-T1)
      ELSE IF (CL.EQ.6) THEN
	   T1=0.5456-0.0236*M1
	   C1=1.-NHCL/8.*(1.-T1)
      ELSE IF (CL.EQ.7) THEN
	   T1=0.8717-0.0179*M1
	   C1=1.-NHCL/8.*(1.-T1)
      ELSE IF (CL.EQ.8) THEN
	   T1=0.9055-0.0638*M1
	   C1=1.-NHCL/8.*(1.-T1)
      ELSE 
	   WRITE(6,*) 'TARKISTA SY™TT™ARVOT!'
	   koodi=1
      ENDIF
      RETURN
      END

      DOUBLE PRECISION FUNCTION TIMECUM(YEAR,iDAY,iMON,iHOUR)
c this function changes the date into the time calculated from the 
c beginning of the year in units "vuorokausi"     
      IMPLICIT DOUBLE PRECISION (A - Z)
      INTEGER YEAR,IDAY,IMON,IHOUR,ID,I
      INTEGER MON(12)
      DATA MON/0,31,28,31,30,31,30,31,31,30,31,30/
      
      ID=IDAY
      if (mod(year,4).eq.0.and.mod(year,100).ne.0.or.
     *mod(year,400).eq.0) MON(3)=29
      DO 172 I=1,IMON
      ID=ID+MON(I)
172   CONTINUE
      TIMECUM=ID+IHOUR/24.                   
      RETURN
      END

      SUBROUTINE DATE(YEAR,TICU,PV,KU,HOUR)
C this function changes the time back to the form day month hour
      IMPLICIT DOUBLE PRECISION (A - Z)
      INTEGER I,YEAR,PV,KU,HOUR
      INTEGER MONTH(12)
      data MONTH/31,28,31,30,31,30,31,31,30,31,30,31/
      
      I=1
      ID=0
      if (mod(year,4).eq.0.and.mod(year,100).ne.0.or.
     *mod(year,400).eq.0) MONTH(2)=29

173   ID=ID+MONTH(I)
      IF (ID.LT.INT(TICU)) THEN
	I=I+1
	GOTO 173
      ENDIF
      PV=INT(TICU)-ID+MONTH(I)
      KU=I
      HOUR=(TICU-INT(TICU))*24+1
      RETURN 
      END

      SUBROUTINE SOLAR(RLAT,RLON,year,iDAY,iMON,HOUR)
C CALCULATION OF SOLAR FLUX S (W/M2), DECLINATION DEK (RAD) AND
C SOLAR HEIGHT ANGLE, GIVEN DATE,TIME (UTC),LAT(deg N),LONG (deg E):
C FORMULAE FROM PALTRIDGE AND PLATT (1977)
      IMPLICIT DOUBLE PRECISION (A-Z)
      INTEGER ALAT,ALON,year,IDAY,IMON,HOUR,ID,I,K,J
      INTEGER MON(12),month(12)
      DATA MON/0,31,28,31,30,31,30,31,31,30,31,30/
      ALAT=NINT(RLAT)
      ALON=NINT(RLON)
      data month/31,28,31,30,31,30,31,31,30,31,30,31/
      
      OPEN(11,FILE='ZENI.DAT', STATUS='UNKNOWN')
      FI=ALAT/57.3
      TICU=TIMECUM(YEAR,IDAY,IMON,HOUR)
      ticu=ticu-2./24.
      call date(YEAR,ticu,iday,imon,hour)
      ID=IDAY-1
      if (mod(year,4).eq.0.and.mod(year,100).ne.0.or.
     *mod(year,400).eq.0) MON(3)=29
      DO 170 I=1,IMON
      ID=ID+MON(I)
170   CONTINUE
      if (mod(year,4).eq.0.and.mod(year,100).ne.0.or.
     *mod(year,400).eq.0) MONTH(2)=29
      DO 300 K=1,month(imon)
      D=ID*6.2832/365.
      DO 200 J=1,24-HOUR
C     S=1368.*(1.00011+.034221*COS(D)+.00128*SIN(D)+.000719*COS(2*D))
      DEK=.006918-.399912*COS(D)+.070257*SIN(D)-.006758*COS(2*D)
     C +.000907*SIN(2*D)-.002697*COS(3*D)+.00148*SIN(3*D)
      H=HOUR+ALON*24/360.
      A=COS(FI)*COS(DEK)*COS((H-12)*6.2832/24)+SIN(DEK)*SIN(FI)
      write(11,107) IDAY,IMON,HOUR,min(90.d0,90.d0-DASIN(A)*57.3)
 107  FORMAT(3I4,1F8.2)
      HOUR=HOUR+1
200   CONTINUE
      iday=iday+1
      HOUR=0
      ID=ID+1
300   CONTINUE
      CLOSE(11)
      END
